ROYAL MCBEE CORPORATION 
ELECTRONIC COMPUTER DEPARTMENT 

I A INTEGRATION OF DIFFERENTIAL EQUATIONS 



(Gill's Method) 
.^ 2 • (Program 28.0) 



FUNCTION I 



To integrate n first-order differential equations of the form 

^= f J (t, y 1 , y 2 , ...,7°). * 

The method is described in detail at the end. (The variables are designated by super- 
scripts because subscripts are reserved for something else.) 

INPUT: 



and m 



Provided in the initial set-up part of the main program, which is used only at 
the very beginning, and in the calling sequence. 

OUTPUT ; 

The new values of the functions, and of other quantities carried along to make 
the calculations possible (see discussion of method). 

For each variable must be stored four quantities: 

i 
y ; 

q\ a. quantity carried along only for use in the computation; 

>f (qfj-qyj) if qf J < Qyj 

I any negative number if q j = q j 

^- 1 y 

q f j > q^J is forbidden. 

q f j is the q of the jth derivative; 

qi is the q of the jth function. 

Note again that the derivative may not be carried at a higher q than the corresponding 
function. 

It is assumed that the functions and derivatives are all carried at the lowest possible 
q's, subject to the restriction of the preceding "paragraph. 

* A later section deals with the problem of integrating equations of order larger 
than 1. 

STORAGE : 

The routine itself occupies 209 cells, (r _l q +0208) 10 ^ but 1 +17-19, L +57, 
L +125, L+153, L +l6l-3 are not used; it also uses 11 cells on tracR 63: 
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6305, 6311, 6318, 6339-1*1, 63U5, 6355-6, 6358, 6362. 

In addition, four quantities must be stored for each variable. 



i variables be stored 


in 


Y 1 


-Y n ; 


the n derivatives 


in 


F 1 


-**,• 


the n "modifiers" 


in 


M 1 


-M*; 


The n q's 


in 


Q 1 


-<f. 



It is seen that, in general, Un storages must be set aside for the functions and the 
quantities necessary to compute them. * 

In the following discussion h is the interval of integration i.e., 
h « t k+1 - t k . 

CALLING SEQUENCE : 

Actually, two calling sequences are necessary. In the "initial set-up" portion of 
the program, used only at the very beginning, must appear ; 

Location Order Address 

<& - 2 B L(C) C is usually B (L +62)_„** 

/* ' 1 H (L o+ lU0) 10 c > L+lh0 10 

^ B L(h) H at (h»l has the hex. repre- 

sentation 7wwwwwwq.) 

^ + 1 H (L o+ 0l5U) 10 > L+15U 

/? + 2 B (L +0160) 10 

/? + 3 C (L o +0038)i > AC 

/? * k H (Q x ) 

/<?* 5 H (Q2) ) q J made for all j. 




A*+ (n+3) H (Q n ) 

* The exception, in integrating n-th- 01 ^ 61 " equations, is discussed later. 
** In one special case this command will be something else. (See page 6.) 

/?and s^+1 store the interval of integration .^+2 and^+3 initialize the routine 
correctly. /9A\ through /# (n+3) make the n fr's 0. T he 6? must be made at the start 
and at any time new initial conditions are read in . The rest of the time they are 
stored by the routine itself, along with the computed values of the yj. The purpose 
of/^+U through ^+( n+3) may, of course, be served by a loop or by data input. 

The independent variable must be carried as the solution of a differential equation. 

% - 1. «(t ) « . . 
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The calling sequence in the main program is: 

Location Order Address 

^ - 15 B L[(n-2) at 29] (n-2) at 29 

H (L o+ oi02) 10 — ► 0102 

B LCM 1 ) M 1 > 

Y (L o+ 0203) 10 A(L o -K)203) 
B LCY 1 ) Y 1 — ► 

Y (L o +0130) 10 A(L o +O130) & 

Y (L o +0131) 10 A(Lo+0131) 
B LCF 1 ) F 1 > 

Y < L o 40 °36) 10 L o+3 6 

B UQ 1 ) Q 1 h 

Y (L o+ 0126) 10 A(L o+0 126) & 

Y (L o+ 0200) 1Q A(L 0+ 200) 

^ - 3 U (c<-l) 

U Special return, for i = h 

R (L o +0l01) 10 

U (L o+ 0200) 10 

<^~ + 1 ordinary return, for i< h 

The return will be too(-2 or to<X+l. Since this is, essentially, a Runge-Kutta 
scheme, there must be 2 different returns. Only values at i = h "count", so the 
"end logic" (testing, printing, etc.) will be done only when i = U. <=<■ +1, therefore, 
will usually be a return to the section of the program in which the derivatives are 
evaluated. 

Notes: It might be thought that h < 1 is a severe restriction; actually, it is not. 
One can easily make a change of variable to get the new interval down to size, and 
change his equations accordingly . 

It is very easy to integrate mth-order equations (m > 1) with this scheme. 
Suppose 

d_y - F(t,y < m - 1 ) , y(m-2) 

— , j } } y) } w here 

dt 

y(j) = d£y_ . 
dtJ 
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Then, y^ 1 ) » fl. 



dt 



y (m) = f 1 = jj (m-D 
dt 

Incidentally, if one stores y, y' 1 ', ... y (m) consecutively, m + 1 cells will suffice 
to hold y and its first m derivatives instead of 2m cells. (Therefore, F^ = y 1 +1.) 
Furthermore, y^)^ 1 < j < m - 1, need be stored only once, not twice. 

ACCURACY : 

Not easy to estimate, but the developer of this method puts the upper bound of 
the error at 

, h y per time step, o r 
_i_ b?ZLy over the whole range. 



TIME: 



Roughly n .1/5 sec# / variable (based on 10.2 msec. / revolution). + h times 
the time to evaluate the derivatives. 

DESCRIPTION OF METHOD ; 

The procedure used is S. Gill's modification of the fourth- order Runge-Kutta 
method, with two changes to enable greater accuracy than provided by equations (26) 
of Gill's paper.* 

For each of n variables must be solved a first- order differential equation 

dyj 



dt 



=fJ to 1 , J 2 , ..., y", t). 



The last equations show what the q J do; it is not yet clear what purposes the modi- 
fiers, mJ = 24 f J " %J 3 serve, h is carried at 0. In (lc), (2c), (3c), and (Uc), 
we have (without superscripts) . 

(5) 7 ± « 7i -l +hr i- 

#S. Gill, "A Process for the Step- by- step Integration of Differential Equations in 
an Automatic Digital Computing Machine." Proceedings Cambridge Philosophical 
Society. Vol Itf, Pt. 1. 
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The r£ are linear combinations of derivatives, and the derivatives and corresponding 
functions may be at different q's. In the machine, then, (5) is calculated by (6): 

(6) 7± at qy = y.^ at q y + (h at ) (1 at ^ . q f) (r A at qf) 

The 1 at (q y -q f ) is 2^-%, or the modifier previously defined. 

Note: It is possible to save (n-1) more cells (and a little time in the integration 
routine) in one special case. Suppose 

^f " % 3 = cons tant f or a H J • 
By far the commonest situation in which this relation would occur is that in which 

q f3 =q Y 3 f ° ra11 *• 

1 
In this case only, only one modifier (in M ) i s necessary, instead of n modifiers 

(in M^M 11 ). Then, one command must be changed and the "preliminary" calling 

sequence must read: 



Location 
^- 2 



Order 

B 
H 

B 



Address 

L(C) In this case, C » U (L q + 1U3) 
(L Q +UiO) 1Q c __^ L+ 1U0 

L(h) 



As used here, the equations for each variable in turn are 



1st step, 






( ^ r l " ^ f o " <tf 
) (lb) «£ = cf 0+ 3rO- 1/2 f Q 

I (1c) y 1 = y + h ri v 



2nd step, 



i = 2 



(2a) rj = (1 -xn72)(f!- q[) 



(2b) 



^2 " 



+ 3r^- (1 - V r T72)(f 1 - q£) 



•1 ""2 

A 



(2c) y 2 = yx+h r 2 



3rd step, 



1 = 3 



(3a) r~ = (1 + VT72T(f 2 - <£) 
(3b) cq = qj + 3rJ - (1 +VI72I f 2 
(3c) y 3 = y 2 + hrj 
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l^th and last step, 



(Ua) r£ = 1/6 f 3 - 1/3 qj 
(Ub) <£-<£ + 3r£ - 1/2 f 3 

(i*c) y^= y 3 + h r J 



One can check out the evaluation of his derivatives by printing the q . The 
following equations are appended for this purpose. 



rjW2 f Q , 

q^f . 



r£~-X>,- 



q^y/T/2 f 



o. 
rj ^ 1/2 f , 

qj ^ 1/2 f Q 

r J /-u j 
q£ r^ . 



If y is linear, the n <~^- " and "^^r" would be "="; qA should always be 
identically 0, and will be except for round-off. 
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STATEMENT OF THE PROBLEM : "Integrate the following system of differential 
equations with the given initial conditions: 

EQUATIONS: ^x = y . dy = _ x . dz = 2 

INITIAL CONDITIONS: x=0 ; y=l ; z=0 

Since the results of the process of numerical integration will give the values 
of the functions or variables, "x","y" &" z " at tne end of a specified time interval, 
an alternative statement of the problem would run as follows: 

"Given the system of differential equations 

dx = „ . dy = -x dz = 2 

oT y > oT ' oT 

and the initial values of the functions at the time t = 0, namely, 

x = ; y = 1 ; z=0 , 

compute the values of the functions at the time t = l/8. 

The punched tape for this problem contains linkage instructions or "calling" 
sequences for three entities: a) functions 

b) derivatives 
I. The problem's "data" c) modifiers 

d) special constants " q " 

II. Subroutine No. 28.0 (Stored on tracks 50, 51, and 52) 
III. Data Output No. 1 (Stored on tracks 08, and 09.) 
The names, initial values, and corresponding scale factors during storage and 
output printing are as follows: 

QUANTITY In. VALUE STORAGE "Q" OUTPUT "Q" LOCATION 

1. x o i h 1003 

2. y '.. 1 1 h 1002 

3. z 7 7 1001 

U. dx/dt 1 h 1035 

5. dy/dt 1 k 103U 

6. dz/dt 2 2 U 1033 

7. $ x 1 1 h 5620 

8. ^ 1 h 5619 

9. q z 2 h 5618 

The operating program may be adjusted, via three instructions, to print the 
above nine quantities at the end of each stage "i", or at the end of a given time 
interval "h", or at the end of a multiple of "h" depending on the words in the 
locations: 2021)., 2027, and 623U. 
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These three printing possibilities may be tabularly listed with their locations and 
corresponding instructions as follows: 

PRINTING AT THE END OF : 

EACH "i" EACH TIME INTERVAL ff h" A MULTIPLE OF "h " 

Location Instruction Location Instruction Location Instruction 

20214. : U 6200 202U : U 6200 202^ : U 1920 

2027 : U 6200 2027 : U 1906 2027 : U 1906 

623U : U 1906 623U : U 1936 623^ : U 1936 

Automatic carriage return key is to be placed such that automatic return of 
the carriage follows the printing of the ninth information item, "q ". (The "q's" 
with the circumflex symbols are special calculational constants and z are not related 
to scaling processes.) 

As can be seen from printing the contents of the problem tape the operating 
program, or control program, is stored in locations on tracks 17, 19, 20, and 62. 
This particular storage arrangement is completely arbitrary. 

A copy of the program description for Subroutine No. 28.0 should be available 
to the reader when using this particular typical problem. 



PROCEDURE 



1. Store Data Output beginning on 0800 

2. Store "A Program for Operating the "Integration of 
Differential Equations" Subroutine No. 28.0"; this 
tape has all of the necessary "Start Fill" Set Mod" 
information all punched on tape. 

3. Store Subroutine No. 28.0 beginning in location 5000. 
h- Execute a transfer to location 1710 via ".0001710' ". 
5. When this is done, the program will compute the 

nine pertinent values and will print them out for 
time equal to l/8 second, as the tape is presently 
punched. To obtain other intervals of printing 
listed above, make changes in words indicated. 



** As the problem tape .is presently punched it will execute 
printing at the end of each time interval "h". 
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Instruction 
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1 


# 












' X 






0.0. °.° 


, ....„,„ ( 


< hk M *Jte*> 




..„> 1 


O.C.Oi o.o o 


/ . . . ! 


1 . . .4- 


' / $.*<? 








2 


. .X/7|£.5.¥-.5 




" 






3 


. . M, 0./ . # o 


' XT. 


; 






4 


, ,X.hl\6.3.4rO 










,0,5 


, , A, -Si 6,3 .(o,i2. 


, a-('i«,<" *M**s*, 








,0,6 


. .KH\£>3.£.2 


T ,L 








,07 


. .A./?, 6.3.6.2. 


<X £ 








8 


. .X.5i6.3.^.<9 


31 *J! 








JL i_ ^ t 9 


. .X./?j £.3.3^ 


/\ (/£/ '« z'A^i tfaS€^ 








,1,0 


, ,//? i^,3 ,5",? , 


£ 








11 


. . .£)£./ .33 


X 








,1,2 


, , .U\O t /.£.(o 










13 


. . ;L/\QO./.+ 










1 4 


, , ,#|C./.^5.7 


Uooz c 








1 S 


, . H\D.O,3.2 


X 








1 6 


, , JJ\OJ ,oj 




9 <3 to 


£*<t 






1 7 












1 8 








, 




,1,9 




X 








2 


. .KM\L.Z.S.L> > 










2 1 


. .X.5 f b.b.U.Z 


4/ 








2 2 


, . .M|aa^.^ 


/WI 7 = c, 








2 3 


. . X ./-/ 1 b.3.5.$ . 


z 




1 ,, 




2 1 


. . «P|0.o. &.o . 


1 

3 









2 5 


, ,X,#,<b.3.4.Z ' 










2 6 


, ,A.Ci 6,3. ^.o - 










2 7 


, , AV^i h,3,5,lo i 


X £ 








2 8 

J>. , 1 1 


. . MJ0.O..5.? - 


/-/■¥ = c, 








L L 2 J 9 

, ,3,0 
. . 3 • l 


_.. ^4>_A^ ' 








, . M\Oj ,3,3 , 


XI 




J 
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* 

Location 


i 

Instruction 

Op. Address 


g- Contents of 
w Address 


Notes 


1 

1 ' 


X! 








1 i i 1 


1 » ' __ 


i i i j 


-1 1 X 


0\O , 3 ( 2 


1 








-*l 1 ,\ 1 ,1- — 


3 3 


. . B\0.i .5.2 


' L^<y(93'? 




■ iii[ 


J 1 ) — . - 


3 4 


. . .N\o.o.3.8 






1 ' i I 




3 5 


1 
. . JJf.i.O. 1 


• X] 


to exit 




•1— — 1 1 , 


3 6 


. . b\l. . .J 






, — i — . — . — f_ 


-1 1 1 


3 7 

- i > 1 — 


. . .U\o.o.z.? 








X.J 1 


3 8 


, . .U^LQO.o.zj 




f~ 1 


: — * * * j 


J 1. . i 


3 9 


. .A A 6.3. i.l 


• K 




■ t i 1 




4 


. .X.S j lo,B .(o.Z 


• 9-z- 




-j « 1 


J 1 1 


4 1 


. . M\0.o.(o.3 


> k^tO+tl) 




- — i — i i (__ 




. 4 2 


. . .Oja/.o.o 


j_ 




'irj 


1 1. -A-— — _ 


4 3 

1 ft .. 1 


. .X.CJ&.3..5:* 


>K 






III, _ 


4 4 

- i i .,t , . 


. .X.5,^J,/./ 


k\ 




1 ■ f 1 


I 1 i ,_., 


4 5 


, . M\0.o.&.3 


i&* i + v~E } 


1 


i , • | 


1 1 I 


4 6 

— ' , ' i... , 


. .X.^J6.3, ij 




1 


■ ■ * 1 i 


1 t 1 - 


4 7 

1 1 i. ., 


. .Xflj 6.3.42 - 


X Si 




— i — i — _j ' 


i 1 i , 


A 8 

.1 — i , I , ... 


, .K4\U3.£2 ' 


— <— a Ji— 




— i — i i 1 i 


... .» i 


4 9 
_ . i i i 


, t XA\Lo,3J.$ . 


-£^^ 




— i — i — i 1 ,i 




5 


. .X.^| 6.3.5.^ . 


^ 




1 L 1 I . | 


• i 


5 1 
_.. i i ■ 


, ,x.fi\k3.s.8 . 


X ^ 


i 


• i . f , 


J,. ...li 


5 2 

« L ( ... 


. . /eJQ/.^.3 - 






1 ' ■ 1 1 


* 1 


5 3 

-i — J~ .„> 


, , .U^.f.^lo . 




( 






5 4 

_ — i — i .. .. 


, , 3\o./.S.c, , 


Uofo3 




1 1 ■ 1 1 


1 — j 


5 5 

: 1... 1 


. . h\g.O.&2 , 


X 




1 1 1 1 ' 


—-j i — 


5 6 


,1 
. . Vf.l.QJ . 




to €.X.tt 


, ■ l—t 1 ' 


— i i 


5 7 

*-- ,1 . , L 


. . . | , , , 






, ; .o.o.Oto. 


°i<>,7 „ 


. , 5 . 8 2 


i 

'!<£" i7. ^-i Sfi b,b ,(o • 


C,*t-y2 




jiil) 


— 1 I , , 


5 9/ 


1 
\^>5,'-5j-S,.5",s5^ ' 


E t 




* i i 1 i 


— t .._ i 


. . 6 .0 A 


..F.rr^F.F.F.F , 


J3 




■ i — — Jt * j,., _ j... 


i 1 _. 


6 1 

— * — -i — . — _ _„ 


i.i!,, fj ' 


3 fc^tf 








6 2 

t__ ...J ._ 1 ._ 

... . 6 . 3 1<£ 


1 

.. .. ... 1 . £. ' 

><■</-,/ 13, j-,j- t A> 


/ 6' £.<j 


M 


- j . .i-_.L. . 


... i l_. i 


_~— zd 1 
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